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ABSTRACT 

In this note we present an approximate method to determine the index of 
refraction of a dielectric obstacle. For simplicity we treat one-dimensional 
models of electromagnetic scattering. The governing equations yield a second 
order boundary value problem, in which the index of refraction appears as a 
functional parameter. The availability of reflection coefficients yield two 
additional boundary conditons. We approximate the index of refraction by a 
k-th order spline which can be written as a linear combination of B-splines . 
For N distinct reflection coefficients, the resulting N boundary value 
problems yield a system of N non-linear equations in N unknowns which are 
the coefficients of the B-splines. 


Research reported in this paper was supported by the National Aeronautics 
and Space Administration under NASA Contract Nos. NAS1-16394, NAS1-17070 and 
NAS1-17130 for the first author and NAS1-17070 for the second author while 
they were in residence at ICASE, NASA Langley Research Center, Hampton, VA 
23665. 



Introduction 


In this note we treat a class of inverse scattering problems in one 
dimension. Similar problems hae been studied and motivated by Hagin [1] and 
Gray and Hagin [2]. Our goal, like theirs is to reconstruct an unknown 
function, which we refer to as the index of refraction, to moderate accuracy 
and by simple means . 

The problem we solve is as follows: 

We seek a function n(x) such that the solution u of the boundary value 
problem 

u"(x) + n^(x)w^u(x) = 0 x€(0,l) 

(u'+in 0 cou) (0) = 2in Q (o (P) 

(u'-in 0 iou) (1) = 0 

satisfies u(0) = f(w) and u'(0) « g(o)), where f and g are known 

functions of the real parameter u). 

This class of problems arises in several situations. In order to 

motivate the above problem we shall describe one such situation arising from 
electromagnetic wave propagation. 

For a general graded index medium, whose refractive index is n(x), the 
time harmonic Maxwell's equation, with single angular frequency k, can be 


written as [3] 
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curl H = ikn^(x)e^E 

curl E = -ikWpH (H) 

o 

div (n (x) g qE) = 0 
div H_ = 0 

where H_ and E_ are the electric and magnetic fields and e ^ and are 

the permittivity and magnetic permeability. Consider a dielectric slab of 
thickness L (Figure 1) parallel to the y-axis 



Figure 1 

We assume the magnetic permeability of the slab is the same as for the ambient 
medium (air). Suppose there is a transverse electromagnetic wave uj 
present. That is to say the electric field _E and the magnetic field are 
at right angles to each other and to the direction of propagation which is 



3 


taken as the x-axis. _E = (0,u(x),0) T satisfies (M)o. H = ( 0, 0, 

satisfies (M^ and (M)^. Demanding II and JE to satisfy (M)^ gives 

2 2 

u"(x) +k p^e^n (x)u(x) = 0. Upon nondimensionalization, in particular 
replacing x with x/L, the above becomes 

u"(x) + (d 2 n 2 (x)u(x) = 0 (1.1) 

2 2 2 

where w = p^e^k L . The bondary conditions are that the tangential 

components of and H_ are continuous across the interfaces of the slab. 

This implies that u and u' are continuous at x « 0 and x « 1. The 

problem that we consider is that the index of refraction of the ambient medium 

is a constant nQ. The process of wave propagation is initiated by a wave 

in^cux 

coming from x = - 00 of the form e • Then the solution of (1.1) for 

x < 0 can be written as 

in n o)x -in^x 

u(x) = e + R(cj)e , (1.2) 

where R(w) is called the reflection coefficient. For x > L part of this 

wave is transmitted and the solution has the form 

in n cox 

u(x) = T(co)e U , (1.3) 

where T(w) is called the transmission coefficient. Since u and u' are 

continuous at x = 0 and x = 1, then (1.2) and (1.3) imply that 

(u'+in^wu) (0) = 2inQ<A) (1.4) 


and 
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(u'-in 0 tou)(l) = 0 (1.5) 

In practice the reflection coefficient R(to) can be measured for any desired 
value of (a). The continuity of u and u' at x = 0 together with (1*2) 
imply that 

u(0) = 1 + R(u)) 

and (1.6) 

u'(0) = in 0 a)(l-R(a))) 

Thus u(0) and u'(0) are known functions of o>. 


2. Solution Procedure 

Our solution procedure consists of approximating n(x) with a k-th order 

spline n(x). This choice of approximation is motivated by the need to solve 

the problem (P) quickly and accurately. We will expound on this statement 

after the description of the scheme. 

We begin with a brief description of the spline spaces S^Cx)* Let 

x - be a partition of [0,1], and let S^(x) denote the space of k-th 

order splines with knots at each x^* s€S k (x) implies that in each interval 

[x ,x ] s is a polynomial of degree at most k-1 and [0, 1] . Let 

r i N*i-V 

{Bf k (x)}_£_^ denote the B-spline basis for S^(x). The general properties 

of B-splines are well known and can be found in [4] for example. A property 

for which we will have occassion to utilize is that ^( x ) = 0 if 

N+K-l 

X ^ X i-k+l ,X i+l^ * ^ a conse( l uence > if fW = L a i B i k (x) 

x€[x ,x ,] then f(x) is simply given by 
m irrri 


and if 
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nrfk-1 


:(x) = l cc^ k (x). 
1=111 9 


( 2 . 1 ) 


The approximation of n proceeds as follows: 
For given N and k we seek an approximation 


n (x) = l X^B 

(n) i=l 1 1,k 


( 2 . 2 ) 


where N = N + k - 1, such that, given the distinct pairs (u^,r(o)j)) 
j = 1, •••,£, where l = f N/2] » the X^^’s are chosen such that the solutions 
u_. of the initial value problems 


u'' + w^n^ (x)u. = 0 
J J (N) J 


u (0) - 1 + R(m ) (I) 


u'(0) = iw (l-R(w )) 


satisfy (u'-iaijU^) (1) = 0 for j = 1,2, •••,£. This is equivalent to solving 

F(X) = 0 (2.3) 

where X. = (x(^,*-«,X (F))T 

and .F ( X.) = ( f i (A) , * * * , f (A) ) T and where 

1 TJ" TT 

f 2j (-) = Real( (uj-iojjTTj) (1)) and = Imag((uj-iw iT ) (1)) 

j = 1,2, •••,£. Before proceeding to the motivation we note that the system 
(2.3) is highly nonlinear. 

In order to solve (2.3) one must solve (I) many times. In general this 
requires the numerical solution of these problems which leads to a large cost 
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in computer time. Our choice of approximation reduces this time to a great 
extent. For example, when k = 1 one can obtain the solution of (I) in 
closed form, and when k > 2 one can easily obtain truncated Taylor 
expansions, whose accuracy is comparable to the numerical integrations of 
(I), but the time taken to compute the expansions is, in general, much less 
than the time taken to integrate (I) numerically. The Taylor expansions are 
computed as follows: For x ^[ x m > x nr fi ] » rT r ,^(x) is a polynomial of degree at 

most k - 1 and is given by 


(*) 


_ m+k-l />— • \ 

n (x) = l Xj B (x). 
(i) i=m 1 ± ’ k 


This follows from (2.1). Let S m (x) and H^x) denote two linearly 

independent Taylor series solutions, expanded about x m , of (I)^. Then for 

xe[x , x . i ] the solution of (I) is given by u. (x) = a S (x) + b H (x) , 

L m m+1 J & J j m m m ™ 


m m 


where the a^s and b m 's are chosen so that u^ and u' are continuous at 
each x m and also to satisfy the initial conditions. 

It is known [4], that if f€(/^[0,l] for j = 0,l, ## *,k-l, then 


inf II f-s II < C, .h j w(f ,h) 

co Ir -i v - 5 J 

ses k (x) K,J 

where h = max (x. ,-x.), C, , is a constant which depends only on k and 
l<i<N 1+1 1 

j, and 


w(f^;^h): = sup{ |f ^ (x)-f ^ (y) | : x,ye[0,l], | x— y | <h} 


is the modulus of continuity of f 


(j) 


estimate cannot be improved. Thus, if n€C v '[0,1], then 


at h. It is also known that the above 
(k) 
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wfn^k ^;^h) < hlln^^ll and the best that we can hope for is lln-n II 

r M 

0[h J. As will be seen, in several numerical examples tested the optimal 
convergence rate was attained. 


3. Numerical Results 


In this section we present some numerical examples for k = 1,2 and 3 
with knot sequence { ^ ^ . In all examples the Levenberg-Marquardt 

algorithm [5] was used to solve (2.3). The r(cj^)'s were computed by 

inputting the exact solution n(x) into (P) • 

To solve a large non-linear system of equations it is usually necessary 

to have a good initial approximation to the solution. For fixed N and k 

our solution was built up according to the following algorithm: 

(i) Select a sequence of integers m = with m^ = k, = N 

and • 

(ii) Compute nr >> by solving (2.3) with an initial guess 


X 




= 1 . 


(iii) For j = 1, 


,M-1 compute n/- \ 

initial guess chosen so that n^ m -j interpolates nj-^ -j 


by solving (2.3) with the 

at the 


m 


j+l J 


distinct points . This gives a linear system for 

's that is invertible if and only if 


the initial X^ ^ 


* - 1 , 


’ m j+i 


[4] . 


x * < ^ < x *+ k 

r ii M 

We found that the best results were obtained by setting m = {2 
for k = 1 and m = { 2i}^ = ^ 


for k = 2 and 3. 
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The maximum attainable N is limited by considerations of computer time 
and accuracy* In the examples tried we found that satisfactory results were 
obtained if N = 32 for k = 1 and N varied between 10 and 14 for k = 
2 and 3 . 

Example 1 * (Figures 2, 6, 10) 
n(x) = 1 + x^ 


Table I 


k 

i 

2 

3 

N 

32 

12 

10 

'"M 1 - 

.036 

.0064 

2.4 x io -6 

computational time (cpu sec*) 

123 

31.2 

8.30 

estimated convergence rate 

1.0 

2.1 



Example 2 * (Figures 3, 7, 11) 

( 1 x€[0,^)u(if, l] 

n(x) =1 

( 1 + sin 2 2ir(x- ■£•) x€[-$>-|] 
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Table II 


k 

i 

2 

3 (a) 

N 

32 

12 

10 

lln-n,. II „ 

(n) 

.16 

.058 

.021 

computational time (cpu sec.) 

195 

101 

35.7 

estimated convergence rate 

1.1 

.96 



Example 3 . (Figures 4, 8, 12) 


n(x) 



xe[o,-j-) 

xe^.i] 


Table III 


k 

i 

2 

3 

N 

16 

8 

10 

«n-n, J\ m 

.50 

.46 

.hi 

M °° 

computational time (cpu sec.) 

25.6 

20.2 

104 


estimated convergence rate 


(a)These results were obtained by demanding the initial approximation to 
interpolate the final approximation for k = 2. 
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Example 4 . (Figures 5, 9, 13) 

{ 1 x€[0,-^) 

■j + 2x xefl , -y] 
2 xe(| ,i] 


Table IV 


k 

i 

2 

3 

N 

32 

14 

10 


.050 

.028 

.025 

computational time (cpu sec.) 

09.2 

117 

51.1 

estimated convergence rate 

1.1 

1.1 



4. Discussion 

In all examples with k = 1 and 2, except for example 3, the optimal 
convergence rate was attained. 

It is obvious from the results listed in Tables 1-4 that if a desired 
accuracy is specified, then the minimum computational time is achieved by 
increasing the order of the splines rather than increasing the number of 
unknowns for a fixed order. 

In practice the measurement of the reflection coefficients is subject to 
experimental error. To simulate this situation we introduced Gaussian type 
noise in the values of R(co). Stable results were obtained if the amplitude 
of the noise was less than 10“^. 
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Our results are comparable in accuracy and computational time to those 
obtained in [1] and [2]* However, our method is applicable to any boundary 
value problem in which the unknown appears as a parameter and sufficient 
additional conditions are known. 
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